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We present a new method of extracting gravitational radiation from three-dimensional numerical 
relativity codes and providing outer boundary conditions. Our approach matches the solution of 
a Cauchy evolution of Einstein's equations to a set of one-dimensional linear wave equations on 
a curved background. We illustrate the mathematical properties of our approach and discuss a 
numerical module we have constructed for this purpose. This module implements the perturbative 
matching approach in connection with a generic three-dimensional numerical relativity simulation. 
QQ ' Tests of its accuracy and second-order convergence are presented with analytic linear wave data. 
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I. INTRODUCTION 



An important goal of numerical relativity is to compute the gravitational waveforms generated by systems of compact 
astrophysical objects such as binary black holes or binary neutron stars. With the prospect that gravitational wave 
detectors such as LIGO, VIRGO and GEO will be on-line in the next few years, it is crucial to study numerical 
relativistic simulations of events which might be observable by these detectors. Such calculations are important not 
only because they could provide signal templates which would considerably increase the probability of detection, but 
also because the comparison of such templates with the observations may provide essential astrophysical information 
on the nature of the emitting sources. The purpose of the Binary Black Hole "Grand Challenge" Alliance 0|, a multi- 
QQ institutional collaboration in the United States, is to study the inspiral coalescence of the most significant source of 
Q\ , signals for the interferometric gravity wave detectors: a binary black hole system. 

' Central to the goal of determining waveforms generated by astrophysical systems is the need for accurate techniques 
which compute asymptotic waveforms from numerical relativity simulations on three-dimensional (3D) spacelike hy- 
I persurfaces with finite extents. In general, the computational domain cannot be extended to the distant wave zone 
^jpj 1^, where the geometric optics approximation is valid. Indeed, computational resource limitations require that the 
• • outer boundary of such simulations lies rather close to the highly dynamical and strong field region, where backscatter 
^ of waves off curvature can be significant. As a result, it is imperative to develop techniques which can "extract" the 
, gravitational waves generated by the simulation and evolve them out to the distant wave zone where they assume 
^ ' their asymptotic form. 

' While the problem of radiation extraction is important for computing observable waveforms from numerical sim- 
ulations, careful implementation of outer boundary conditions is also crucial for maintaining the integrity of the 
simulations themselves, as poorly implemented boundary conditions are a likely source of numerical instabilities. 
These outer boundary conditions are also decisive in framing the desired physical context for the simulation, e.g., an 
isolated source in an asymptotically fiat spacetime. For typical applications, we can summarize the requirements of a 
radiation-extraction/outer-boundary module as: (a) supporting stable evolution of Einstein's equations, (b) minimiz- 
ing spurious (numerical) reflection of radiation at the boundary, (c) providing accurate and numerically convergent 
approximations to the gravitational waveforms that would be observed in the wave zone surrounding an isolated 
source, (d) incorporating effects of radiation reflection off background curvature outside the numerical boundary 
when appropriate (for example when the outer boundary is in a strong fleld region). 

In this paper we present a new method for extracting gravitational waveforms from a 3D numerical relativity code 
while simultaneously imposing outer boundary conditions. Our approach is motivated by earlier investigations of 
gauge-invariant extraction techniques , but promises to be more generally applicable in cases where the background 
curvature is significant near the outer boundary of the computational 3D grid. Our method matches a full 3D Cauchy 
solution of Einstein's equations on spacelike hypersurfaces with a perturbative onc-dimensional (ID) solution in a 
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region where the waveforms can be treated as hnear perturbations on a spherically symmetric curved background 

The plan of this paper is as follows: in Section ||, we describe the mathematical basis of our method and derive 
the linearized radial wave equatio ns w hich account for the evolution of the gravitational waves in the perturbative 



region of the spacetime. In Section III , we discuss the strategies for the numerical solution of the above equations and 
present a numerical code we have constructed which represents a general module implementing our extraction/outer 
boundary method in conjunction with a 3D numerical relativity simulation. In this paper we focus on tests of the 
outer boundary module as a self-contained unit using analytic solutions. In companion papers we focus on tests of 
the module in the practical context of a typical application. In |^ , we have presented tests of this outer boundary 
module in conjunction with the 3D "interior code" of the Alliance which evolves Einstein equations in the standard 
"3 + 1" form (as presented in ||^). A more thorough discussion of these results is forthcoming 0. 

II. THE CAUCHY-PERTURBATIVE MATCHING METHOD 

Einstein's equations are highly nonlinear and when spacetime is characterized by rapidly varying strong fields, 
the full 3D nonlinear equations must be used. Outside of an isolated region of this kind, however, a perturbative 
approximation, in which gravitational data are treated as linear perturbations of an exact solution to Einstein's 
equations, may be valid. In this perturbative region a linearized approximation to Einstein's equations could then 
exploited to simplify the evolution of gravitational data. 

The idea behind a Cauchy-perturbative matching approach is to supplement the computationally expensive evolu- 
tion of the full Einstein equations with the comparatively simpler evolution of the linearized equations in a perturbative 
region. Figure |l] provides a schematic picture of the Cauchy-perturbative matching approach. The square region cov- 
ered by the grid represents the 3D computational domain N (one dimension is suppressed) on which Cauchy evolution 
of the full Einstein equations is computed. The dark central area in N includes the strong field highly dynamical region, 
where the nonlinear Einstein equations must be solved. The medium and light shaded annular area, T', represents the 
perturbative region. Anywhere in the (medium shaded) intersection of N and V, we can place an extraction 2-sphere 
E, of radius r^, where the gravitational field information is read out. This information is then evolved (by means of 
the linearized Einstein equations) in V, which ranges from E out to a large distance (shown as a dotted circle) where 
the asymptotic waveforms can be identified. Outer boundary data for N can be constructed from perturbative data 
in the intersection of N and V. 

Previous investigations j^] achieved the desired perturbative simplification by matching the nonlinear solution onto 
analytic solutions of Einstein's equations linearized on a Minkowski background. Further simplification was achieved 
by decomposing perturbative data in a multipole expansion. We extend this approach to cases where curvature is 
significant by choosing as our approximation a linearization of Einstein's equations on a Schwarzschild background. 
In principal one could generalize further to a Kerr background, and this will be the subject of future work. We also 
decompose perturbative data on this background with a multipole expansion, reducing the 3D linearized equations 
to a set of ID equations for each multipole mode. This reduction allows us to evolve data everywhere in P on a 
one-dimensional grid, L. It is important to note that all of the 3-dimensional tensor data in V can be reconstructed 
from the multipole amplitudes on L. Our method, therefore, is to match a computationally expensive evolution of 
the full Einstein equations onto a considerably less expensive evolution on a ID grid in a region where background 
curvature is still significant. 



'^An alternative approach to the problem of wave extraction and outer boundary conditions has been developed to match the 
Cauchy solution to solutions on characteristic hypersurfaces 
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FIG. 1. Schematic of matching procedure for two successive timeslices (one dimension is suppressed). N is the 3D compu- 
tational domain. The dark shaded region shows the strong field highly dynamical region in N. The medium and light shaded 
annular region represents the "perturbative region", V, where the linearized equations are evolved on a ID grid L (not shown). 



A. Hyperbolic formulation 

Rather than characterize radiation asymptotically in terms of certain variables constructed from the metric |^ , we 
use a new approach which characterizes radiation in terms of the extrinsic curvature. This is made possible by a 
recently developed spatially gauge-covariant hyperbolic formulation of general relativity. This system is constructed 
from first derivatives of the spacetime Ricci tensor and may therefore appropriately be called the "Einstein- 

Ricci" system. 

The Einstein-Ricci equations are obtained from the "3 + 1" form of the metric, 

ds^ = -N^dt^ + g,j{dx' + I3'dt){dx' + (3^ dt) , (1) 

where N is the lapse function, /3* is the shift vector, and gij is the spatial metric in the slice S. An appropriate time 
derivative operator that evolves spatial quantities along the normal to the slice S is 

do^dt~Cp , (2) 

where Cp is the Lie derivative along the shift vector in E. 
The extrinsic curvature Kij of S can be defined by 

^og^J = -2NK,j , (3) 

which serves also as the evolution equation for the spatial metric. By working out the expression = 9o-Rij — 2 V(ji?j)o 
in 3 + 1 form, where Rij and i?jo are components of the spacetime Ricci tensor and Vi denotes the spatial covariant 
derivative, we find a wave-like equation which governs the evolution of Kij: 

- NaKij = Jy + Sij - , (4) 
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where the physical wave operator for arbitrary shift is □ = —N~^doN~^do + VfcV'^. Equation (Q) is an identity until 
we substitute the Einstein equations Raf3 = STrlTa/s — ^T^gap) into ^iji G = c = 1 ). 

The detailed form of the right hand side of (||) can be found in |^,|lO| ; the present conventions arc those in [0 . 
Here we simply point out that Sly has become a matter source that is zero here, Jij is the nonlinear self-interaction 
term in 3 + 1 form, and Sij is a slicing-dependent term that must involve fewer than second derivatives of Kij to 
render a true (hyperbolic) wave equation. A simple way to satisfy the restriction on Sij is to invoke the harmonic 
slicing condition 

doN + N'^H = 0, (5) 

where H is the trace of -ft'ij, and from which follows Sij — 0. 

For appropriate choice of initial data , equations (||) , , and (|^) represent the dynamical part of Einstein's 
equations. Combining them wc obtain a quasi-diagonal hyperbolic equation for f;^, with principal (highest-order) 

part uBq. Hence (|), (|), and (|) may be said to give the "third-order" form of the Einstein-Ricci system. 

We note that the third-order Einstein-Ricci system can also be cast into a first-order symmetric hyperbolic form 
[p[-]ri|. 1^ It also possesses a higher order form (the "fourth-order Einstein-Ricci system"), essentially a wave equation 
for {doKij), obtained from doilij + ViVji?oo This system has a well-posed Cauchy problem and complete 

freedom in choosing both /3* and N: it has no analog of a slicing term like Sij. This fourth-order form is used to 
develop fully gauge-invariant perturbation theory in 

B. Perturbative Expansion 

The first step in obtaining radial wave equations is to linearize the hyperbolic Einstein-Ricci equations around a 
static Schwarzschild background. We separate the gravitational quantities of interest into background (denoted by a 
tilde) and perturbed parts: the 3-metric gij — gij + hij, the extrinsic curvature Kij — Kij +Kij, the lapse N = N + a, 
and the shift vector — + u*. In Schwarzschild coordinates {t, r, 0, ip), the background quantities are given by 



r ) 

g.jdx'dx^ = N-^dr^ + r^{d0^ + sin^ Odcj)^) , (6b) 

f = , (6c) 

i^y = , (6d) 

while the perturbed quantities have arbitrary angular dependence. The background quantities satisfy the dynamical 
equations dtgij = 0, dtN = 0, and thus remain constant for all time. The perturbed quantities, on the other hand, 
obey the following evolution equations 

dthij = -2NKij + 2V(jUj) , (7a) 
dta = v'V.N - N^K , (7b) 
N~^dtK^j - iVV'^VfcK.j = - 4V(jk''^-) VfcTV -I- N^^K^jV^NVkN + iV'^NVkHtj 

+ 4:d(^,Kdj)N + 2N-^KV,N^jN - 2iV^fc(,K^) - 2NRkrj„rK^'^ , (7c) 

where k = and the tilde denotes a spatial quantity defined in terms of the background metric, Ijij. Note that the 
wave equation for Kij involves only the background lapse and curvature. 



^In the equation for doT^jk was inadvertently omitted. See 
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C. Angular decomposition 



We can further simplify the evolution equation ^idj) by separating out the angular dependence, thus reducing it to 
a set of ID equations. We accomplish this by expanding the extrinsic curvature in Regge- Wheeler tensor spherical 
harmonics Jl6t and substituting this expansion into (^) . Using the notation of Moncrief we express the expansion 
as 



Kij = a^{t,r){ei)ij + rb^{t,r){e2)ij + N ^a+(i, r)(/2).y + r6+(i, r)(/i),y + 



(8) 



where (ei)^ , • • • , (/4)ij are the Regge- Wheeler harmonics, which are functions of [9, (p) and have suppressed angular 
indices {i,m) for each mode. The odd-parity multipoles (ax and bx) and the even-parity multipoles (a+, 6+, c+, and 
also have suppressed indices for each angular mode and there is an implicit sum over all modes in (H). The six 
multipole amplitudes correspond to the six components of Kij. However, using the linearized momentum constraints 



(9) 



we reduce the number of independent components of to three. An important relation is also obtained through 
the wave equation for k, whose multipole expansion is simply given by k = h{t, r)Y^^ where Y^^ {9, 0) is the standard 
scalar spherical harmonic and again there is an implicit sum over suppressed indices {i,m). Using this expansion, 
in conjunction with the momentum constraints we derive a set of radial constraint equations which relate the 
dependent amplitudes (&x)f,„, {b+)t^, and to the three independent amplitudes (ax)f„, (a+)<.„, , (^)fm- 
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(10a) 
(10b) 
(10c) 
(lOd) 



for each (I, to) mode. 

Substituting (^) into (^) and using the constraint equations ( p^ ) , we obtain a set of linearized radial wave equations 
for each independent amplitude. For each (i?, to) mode we have one odd-parity equation 
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and two coupled even-parity equations. 
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These equations are related to the standard Regge- Wheeler and Zcrilli equations ||l6|,|l7|], which can be derived in a 
more complete analysis of gauge- invariant hyperbolic formulations JT^ . 

The radial wave equations ([ll|)-(|l3|) for each {£, m) mode of the independent multipole amplitudes (ax);„, (a-i-)f„, 
(/i)^^ form the basis for our approach. In the perturbative region, they replace the nonlinear Einstein equations and 
determine the evolution of Kij. They can be used to evolve, with minimal computational cost, gravitational wave 
data to arbitrarily large distances from the highly dynamical strong field region. The evolution equations for hij ( [7a| ) 
and a (7b) can also be integrated using the data for Kij computed in this region. Note that because hij and a evolve 
along the coordinate time axis, these equations need only be integrated in the region in which their values are desired, 
not over the whole region L (these quantities have characteristic speed zero). 
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III. NUMERICAL IMPLEMENTATION 



This section is a general guide for the numerical implementation of the Cauchy-perturbative matching method for 
radiation extraction and outer boundary conditions described so far. 

Consider a 3D numerical relativity code which solves the Cauchy problem of Einstein's equations in either the 
standard ADM form or in the hyperbolic form JT^ . During each timestep the procedure followed by our module 
for extracting radiation and imposing outer boundary conditions can be summarized in three successive steps: (1 ) 
extraction of the independent multipole amplitudes on E, (2) evolution of the radial wave equations (|l|)-(|l|) on L 
out to the distant wave zone, (3) reconstruction of Kij and dtKij at specified gridpoints at the outer boundary of N. 
We now discuss in detail each of these steps: 

(1 ) Extraction 

As mentioned in Sect. the extraction 2-sphere E acts as thejoining surface between the evolution of the 
highly dynamical, strong field region (dark shaded area of Fig. and the perturbative regions (light shaded 
areas). At each timestep, Kij and dtKij are computed on N as a solution to Einstein's equations. In the 
test cases presented here, N uses topologically Cartesian coordinates, although there are no restrictions on the 
choice of the coordinate system. The Cartesian components of these tensors are then transformed into their 
equivalents in a spherical coordinate basis and their traces are computed using the inverse background metric, 
i.e. H = g^^Kij, dtH = g^^dtKij. From the spherical components of Kij and dtKij, the independent multipole 
amplitudes for each (£, m) mode are then derived by an integration over the 2-sphere: 

(«x = I [Kr^ do ~ Kre d^] dn , (14a) 

= / dn , (14b) 

(/i),„ = [ H Yl^dn . (14c) 



Their time derivatives are computed similarly. Rather than performing the integrations (14a)-(14c) using 
spherical polar coordinates, it is useful to cover E with two stereographic coordinate "patches". These are 
uniformly spaced two-dimensional (2D) grids onto which the values of Kij and dtKij are interpolated using 
either a three-linear or a three-cubic polynomial interpolation scheme. Each point on the 2-sphere, denoted by 
spherical coordinate values (0,0), corresponds to a point {q^p) on a stereographic grid whose coordinates can 
be combined into a single complex number <^: 



tan I 2 ' ^ 



i0 



(15a) 
(15b) 



where N and S denote the northern (0 < < 7r/2) and southern {n/l <9 <n) hemispheres, respectively. As a 
result of this transformation, the integrals over the 2-sphere in (|lj) are computed over the stereographic patches, 
which naturally avoid polar singularities (see |20) for a complete discussion of the properties and advantages of 
the stereographic coordinates). In our tests, the integrals over each patch are computed using a second-order 
stereographic quadrature routine developed within the Alliance pO| . 

(2) Evolution 

Once the multipole amplitudes, (ax)f„) (o-i-)f„, (^){„ and their time derivatives are computed on E in the 
timeslice t = to, they are imposed as inner boundary conditions on the ID grid. Using a second-order integration 
scheme (we have tested both Leapfrog and Lax-Wendroff ^^), our module then evolves the radial wave equations 
(|ri|)-(P^ for each {£,m) mode forward to the next timeslice at i = ti. The outer boundary of the ID grid is 
always placed at a distance large enough that background field and near-zone effects are unimportant, and a 
radial Sommerfeld condition for the wave equations (pT|)-(p^ can be imposed there. Of course, the initial data 
on L must be consistent with the initial data on N. This can either be imposed analytically or determined by 
applying the aforementioned extraction procedure to the initial data set at each gridpoint of L in the region of 
overlap with N. In the latter case, initial data outside the overlap region can be set by considering the asymptotic 
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fall-ofF of each variable. It should be noted that in the Cauchy-characteristic matching approach initial data also 
must be set in the characteristic hypersurfaces and, for realistic sources like binary black holes, will necessarily 
be approximate. 

(3) Reconstruction and Matching 

From the perturbative data evolved to time ti, outer boundary values for N can now be computed. The procedure 
for doing this differs depending on whether a hyperbolic or an ADM formulation of Einstein's equations is used 
by the 3D "interior code". For a hyperbolic code (cf. it is necessary to provide boundary data for Kij and 

dtKij. For an ADM code (cf. on the other hand, outer boundary data only for Kij are necessary, since 

the interior code can calculate gij at the outer boundary by integrating in time the boundary values for Kij . In 
either case, if outer boundary values for the lapse N are needed [e.g. for integrating harmonic slicing condition 
], these can be computed by the perturbative module or by integration of H at the boundary. 

In order to compute Kij at an outer boundary point of N (or any other point in the overlap between N and 
V 01), it is necessary to reconstruct Kij from the multipole amplitudes and tensor spherical harmonics. The 
Schwarzschild coordinate values {r,9,(j)) of the relevant gridpoint are first determined. Next, (ax)(,„, (0+)^^, 
and (h^^) for each {i,m) mode are interpolated to the radial coordinate value of that point. The dependent 
multipole amplitudes {bx)i^, and ((1+)^,^ are then computed using the constraint equations (px|). 

Finally, the Regge- Wheeler tensor spherical harmonics (ei)ij-(/4)ij are computed for the angular coordinates 
{9, 4>) for each (^, to) mode and the sum in equation (^ is performed. This leads to the reconstructed component 
of Kij (and therefore Kij)\ a completely analogous algorithm is used to reconstruct dtKij. 

It is important to emphasize that this procedure allows one to compute Kij at any point of N which is covered 
by the perturbative region. As a result, the numerical module can reconstruct the values of Kij and dtKij on 
a 2-surface of arbitrary shape, or any collection of points outside of E. 

Numerical implementation of this method is rather straightforward. Very few modifications to a standard 3D 
numerical relativity code are necessary in order to allow for the simultaneous evolution of the highly dynamical 
region and of the perturbative one. Because of the use of numerically inexpensive integration of ID wave equations, 
implementation of this module provides gravitational wave extraction and stable outer boundary conditions with only 
minimal additional computational cost. 

Finally, it should be noted that, in practice, we may not know a priori if the Schwarzschild-perturbative approxima- 
tion is valid near the outer boundary of a given numerical relativity simulation. Through experimentation, however, 
it is possible to test the validity of the approximation. This can be done, for instance, by extracting data at different 
radii and comparing the waveforms computed at the outer sphere with those evolved from the inner sphere. This 
makes it possible to determine if the neglected terms in the approximation have a significant effect. At any point in 
the overlap region between N and V, it is possible to reconstruct gravitational wave data and compare these values 
with those computed by the full nonlinear evolution. 

IV. NUMERICAL TESTS 

In order to establish the accuracy and convergence properties of our code we have studied the propagation of linear 
waves on a Minkowski background (M = 0) . This is a natural first test since we can compare each stage of the 
numerical procedure described in Section |l| against a known analytic solution |^,^ . 

In these tests we assign analytic values to each gridpoint of N at every time step. This allows us to study the 
accuracy and convergence properties of the module independently of any errors which may develop in a 3D numerical 
evolution of linear waves. Elsewhere 0, we will present results of tests of this module running with a full 3D evolution 
code (i.e. the interior code of the Alliance ||l^), with emphasis on the issues of stability of the outer boundary and 
accuracy of extracted waveforms. 

We have considered analytic data for i — 2, m — Q even-parity linear waves, initially modulated by a Gaussian 
envelope with amplitude A = 10~^ and width parameter 6=1. These waves are time-symmetric at t = and 
thus have ingoing and outgoing parts. The 3D grid is vertex-centered with extents {x,y,z) £ [—4,4] and resolutions 
ranging from (17)^ to (129)"^ points [corresponding to (16)^ and (128)'^ zones, respectively]. The resolution of the 
stereographic coordinate patches corresponds to the resolution of N and therefore ranges from (16)^ to (128)^ zones 
on each hemisphere. For the specific tests presented here, E is located at a radius r^ = 3 and similar results have 
been obtained also for r^ = 0.5,1.0,1.5,2,2.5,3.5. In fact, on a fiat background spacetime and for weak waves on 
Schwarzschild-like backgrounds, the perturbative approximation is valid throughout the 3D domain and the position 
of E is thus arbitrary. 
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Since these waves are traceless and even-parity, with pure £ = 2, m = angular dependence, the only non-zero 
independent multipole amplitude we expect to find at E is (a+)2Q. Diagram (a) of Fig. ^ shows plots of {a+)^g 
extracted at function oi t — r for various resolutions of N. The amplitude is scaled by to compensate 

for the radial fall-off. 

The curves in diagram (a) clearly show that the extracted waveform approaches the analytic value (denoted by 
a solid line) as the resolution is increased. However, in order to establish the exact rate at which the computed 
solution approaches the analytic one, we have also performed convergence tests. These tests are designed to check 
that no coding error has been made and that the numerical scheme employed in the solution is providing results at 
the expected accuracy. While there are a number of different ways to perform these tests, we have exploited the 
knowledge of an analytic solution and computed the residuals R between the computed solution and the analytic 
one J-'a as a function of the resolution (or, equivalently, of the number of gridpoints). For a second-order accurate 
numerical scheme (as the one used here) on a uniform cubical grid, we expect the residuals to follow the simple law 

RiN^) ^T, -Ta = 0{h^) , (16) 

where 0(Y?) contains the second and higher order error terms and h — L/{N — 1) is the grid resolution, with L 
being the spatial dimension of the grid. If the numerical computation is second-order accurate and a number of 
simulations with different grid resolutions, each differing by a factor 2, are performed, we should expect the residual 
to fall quadratically to zero. Diagram (b) of Fig. ^ shows this is indeed the case; there, we have multiplied the 
residuals obtained with different resolutions by the coefficients that make the leading order error terms comparable. 
The good overlapping of the different curves is an indication that a second-order convergence has been achieved. 




FIG. 2. (a) Timeseries of the multipole amplitude (a+)2o extracted at a 2-sphere of radius = 3 for various grid resolutions. 
The amplitude is scaled by r'^ to compensate for the radial fall-off. (b) Residuals of the leading order error term for different 
grid resolutions differing by a factor of 2. The residuals are multiplied by 4, 16, and 64 in order to make the errors comparable. 
If no higher-order terms were present, all of the curves would coincide. 

The accuracy of this extraction procedure can also be tested by examining the waveforms for the other multipole 
amplitudes computed which analytically vanish. Figure |^ shows plots of several even-parity (upper diagram) and 
odd-parity (lower diagram) amplitudes computed at the extraction 2-sphere for an resolution in N of (65)'' points. 
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As a result of numerical truncation error introduced in the extraction procedure, these modes are not exactly zero. 
However, even the largest amplitude mode is over three orders of magnitude smaller than the only analytically non- 
vanishing independent amplitude (a+)2o. Moreover, all of these amplitudes are second-order convergent to zero as 
the resolution is increased. Similar considerations apply also for the (h)^^^ multipole amplitudes: although the data 
is analytically traceless, very small (h)^^ modes are extracted at the 2-sphere. These modes, which we will not show 
here, are the order of round-off error (approximately 10^^^ for these tests) and may be considered as effectively zero. 




012345678 
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FIG. 3. Timeseries of the analytically vanishing even (a+)j^ and odd (ax)f„ parity multipole amplitudes. The extraction 
is made at a 2-sphere of radius — S and N has (65)^ points. 

Next, we consider the accuracy of the evolution in the perturbative region of the extracted amplitudes. The time 
integration of (|Tl|)-([l^) on L is performed using a Leapfrog integration scheme with a spatial resolution adjusted so 
that the timesteps in N and L are identical. This imposes a relation involving the gridspacing of N and the ratio of 
Courant factors for N and L. Such a choice ensures a correspondence between resolutions of N and V Fig. ^ shows 
plots of {a-^-).2g evolved to a radius r = 8 from the extracted signal at r = = 3. Different curves correspond to 
different resolutions and show the convergence to the analytic solution. The outer boundary of L is located at r = 33, 
where outgoing wave Sommerfeld conditions are imposed. For radial scalar wave equations, this represents a very 
good approximation which has been shown to be both accurate and stable. 
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FIG. 4. Timeseries of (a+)2o evolved to r = 8 for various grid resolutions. Here also, the amplitude is scaled by to 
compensate for the radial fall-off. 

Finally, we consider the accuracy in the reconstruction of the outer boundary data. Since we are using analytic 
data in N, we can only compare the outer boundary data with the analytic ones. In a forthcoming paper 0], where 
we will make use of a numerical solution of Einstein's equations, we will also discuss the issues of stability related to 
the use of a Cauchy-perturbative matching method. For conciseness we consider here reconstructed outer boundary 
data only for Kij; the reconstruction of dtKij follows analogously. Diagram (a) of Fig. ^ shows the timeseries of the 
reconstructed value of K^z computed at the point (a; = 4, y = 0, z = 0) for various resolutions and its comparison 
with the analytic solution. Also in this case, diagram (b) of Fig. |^ gives proof of the second-order convergence of 
the numerical module even if, in this case, higher order error terms become apparent with the very coarse resolution 
simulations [i.e. in the case of (17)^ gridpoints]. The small peak observed at t — r ss 1 is the result of a slight difference 
between the analytic initial data on L and the extracted signal at t = 0. This error rapidly disappears as the resolution 
is increased. 
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FIG. 5. (a) Timeseries of the reconstructed values for K^z at the grid point (4, 0, 0) for various grid resolutions, (b) 
Residuals of the leading order error term for different grid resolutions differing by a factor of 2 (cf. Fig. 0). 

A more global measure of the accuracy and of the convergence properties of the boundary data is obtained by 
computing the L2 norm of the error in Kij as measured over the whole 3D outer boundary. In Fig. |^ we plot the 
L2 norms of K^z at successive resolutions, normalizing these differences by the factor which would make the plots 
overlap if the convergence to analytic data were exactly second-order. Here we again see that the desired convergence 
rate is achieved over the whole boundary, particularly at finer grid resolutions. 
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FIG. 6. Plot of L2 norms of K^z computed over the outer boundary for successive grid resolutions. The norms at different 
grid resolutions are scaled so that they overlap if truly second-order convergent. 



V. CONCLUSION 

We have presented a method for matching gravitational data computed from a 3D Cauchy evolution of Einstein's 
equations to a computationally simpler evolution of radial wave equations linearized on a Schwarzschild background. 
This method should be applicable to a variety of physical problems where curvature is significant throughout the 
computational domain, so long as the time-dependent fields can be treated as linear perturbations on a spheri- 
cal background. Our approach promises to offer an accurate means of computing asymptotically gauge-invariant 
waveforms at large distances from the domain of the simulation and to provide stable, physically correct boundary 
conditions. 

We have also discussed a numerical code we have developed that implements this procedure and can be used with 
a general numerical 3D simulation, solving Einstein's equations in either the hyperbolic or ADM formulation. This 
code correctly extracts waveforms from analytic linear wave data and recomputes that data at the outer boundary of 
the 3D grid. A more extensive discussion of the stability properties of this approach will be discussed in a forthcoming 
paper |Q, as well as practical issues arising from application to a real evolution code environment. 

Our Cauchy-perturbative matching method can be extended to more general circumstances, e.g. perturbations on 
axially symmetric backgrounds or other slicings of Schwarzschild black holes. Similar analyses using other hyperbolic 
formulations [l^j2^ may also provide important insight to the physical understanding of radiation extraction and 
lead to modules which work with simulations based on these formulations. 
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